#ifndef FEMSPACE_H
#define FEMSPACE_H


#include "Mesh.h"
#include "Element.h"
#include "PoissonEquation.h"
#include "BoundaryFunction.h"
#include <map>
#include <Eigen/Sparse>
#include "Eigen/IterativeLinearSolvers"
typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::Triplet<double> Tri;

/**
 * @brief 有限元空间基类，各种当方程的求解可以继承该类。
 *        提供的功能有：组装刚度矩阵，组装右端项，边界条件处理，计算数值解并返回；
 * 
 */
template<unsigned int DIM>
class FEMSpace
{
public:
    /** 
     * 默认构造函数
     * 
     */
    FEMSpace() = default;
    /**
     * @brief 组装刚度矩阵；
     */
    virtual void AssembleStiffMatrix() = 0;
    /**
     * @brief 组装右端项；
     * 
     */
    virtual void AssembleRightHandsTerm() = 0;
    /**
     * @brief 处理边界条件
     * 
     */
    virtual void DealWithBoubdaryCondition() = 0;
    /**
     * @brief 求解
     * 
     */
    virtual void Solve() = 0;
    /**
     * @brief 返回解
     * 
     * @return Eigen::VectorXd 
     */
    virtual Eigen::VectorXd solution() = 0;
    /**
     * @brief 局部刚度矩阵在i,j处的值
     * 
     * @param i 下标１
     * @param j 下标２
     * @return double 值
     */
    virtual double LocalStiffMatrix(int i, int j) = 0;
    /** 
     * 计算L2误差
     * 
     * @param real_solution 准确解
     * 
     * @return L2误差
     */
    virtual double Error_L2Norm(std::function<double(double*)> real_solution) = 0;
     /** 
     * 计算H1误差
     * 
     * @param real_solution 准确解
     * 
     * @return H1误差
     */
    virtual double Error_H1Norm(std::function<double(double *)> real_solution,std::function<double(double*)> real_solution_x,std::function<double(double *)> real_solution_y) = 0;
    
    /** 
     * 默认析构函数
     *  
     */
    virtual   ~FEMSpace() = default;
};


class Poisson_2D:public FEMSpace<2>
{


public:
    Poisson_2D();
    /** 
     * 指定网格、单元、方程、边值条件
     * 
     * @param mesh 网格
     * @param element 单元
     * @param equation 要解决的方程
     * @param BoundaryCondition 边值条件 
     * 
     */
    Poisson_2D(Mesh<2>* mesh,Element<2>* element,PoissonEquation* equation,BoundaryCondition<2>&BoundaryCondition);
    /** 
     * 拷贝构造函数
     * 
     * @param p Poisson_2D资料
     * 
     */
    Poisson_2D(const Poisson_2D& p);
    ~Poisson_2D() = default;
    virtual void AssembleStiffMatrix();
    virtual void AssembleRightHandsTerm();
    virtual void DealWithBoubdaryCondition();
    virtual void Solve();
    virtual Eigen::VectorXd solution();
    virtual double LocalStiffMatrix(int i, int j);
    virtual double Error_L2Norm(std::function<double(double*)> real_solution);
    virtual double Error_H1Norm(std::function<double(double *)> real_solution,std::function<double(double*)> real_solution_x,std::function<double(double *)> real_solution_y);
    /** 
     *返回刚度矩阵
     * 
     * 
     * @return Eigen::SparseMatrix<double>
     */
    SpMat A() const;
    /** 
     * 返回右端项
     * 
     * 
     * @return Eigen::VectorXd
     */
    Eigen::VectorXd Rhs() const;
protected:
    SpMat _A;
    Eigen::VectorXd _Rhs;
    Eigen::VectorXd _Solution;
    Mesh<2>* _mesh;
    Element<2>* _element;
    PoissonEquation*_equation;
    BoundaryCondition<2> _boundarycondition;
    
};
Poisson_2D::Poisson_2D()
{   
    _A = SpMat(int(_mesh->n_dofs()), int(_mesh->n_dofs()));
    _Rhs = Eigen::VectorXd(int(_mesh->n_dofs()));
}
Poisson_2D::Poisson_2D(Mesh<2>* mesh,Element<2>* element,PoissonEquation* equation,BoundaryCondition<2>& BoundaryCondition)
{
    _mesh = mesh;
    _element = element;
    _equation = equation;
    _boundarycondition = BoundaryCondition;
    _A = SpMat(int(_mesh->n_dofs()), int(_mesh->n_dofs()));
    _Rhs = Eigen::VectorXd(int(_mesh->n_dofs()));
}

Poisson_2D::Poisson_2D(const Poisson_2D& p)
{
    _mesh = p._mesh;
    _element = p._element;
    _equation = p._equation;
    _boundarycondition = p._boundarycondition;
    _A = p._A;
    _Rhs = p._Rhs;
}

void Poisson_2D::AssembleStiffMatrix()
{
    int n_Dofs = _element->n_Dofs();
    std::vector<Tri> TriList((_mesh->n_element()) * n_Dofs * n_Dofs);
    std::vector<Tri>::iterator it = TriList.begin(); 
    /// traverse each element
    for (int k = 0; k < _mesh -> n_element(); k++)
    {
        std::vector<int> Idx = _mesh->NodeofEle(k);
        std::vector<Dofs<2> > temp(n_Dofs);
        for(int i = 0;i < n_Dofs; i++)
        {
	    //temp[i].GlobalIndex() = Idx[i];
	    temp[i] = _mesh->DofsofIndex(Idx[i]);
        }
        _element->SetDofsList(temp);
        for (int i = 1; i <= n_Dofs; i++)
	    for (int j = 1; j <= n_Dofs; j++)
	    {
		*it = Tri(_element->NdIdx(i), _element->NdIdx(j), LocalStiffMatrix(i, j));
		it++;
	    }
    }
     _A.setFromTriplets(TriList.begin(), TriList.end());
     _A.makeCompressed();
     //std::cout << _A << std::endl;
}

void Poisson_2D::AssembleRightHandsTerm()
{
    int n_Dofs = _element->n_Dofs();
    int n_GaussPnt = _element->n_GaussPnt();
    auto _func  = _equation->RightHandsTermFunction();
    for (int k = 0; k < _mesh -> n_element(); k++)
    {
        std::vector<int> Idx = _mesh->NodeofEle(k);
        std::vector<Dofs<2>> temp(n_Dofs);
        for(int i = 0; i < n_Dofs; i++)
        {
	    //temp[i].GlobalIndex() = Idx[i];
	    temp[i] = _mesh ->DofsofIndex(Idx[i]);
        }
        _element->SetDofsList(temp);
        for(int i = 1; i <= n_Dofs; i++)
        {
            double a = 0.0;
            for (int j = 0; j < n_GaussPnt; j++)
            {
		        double xi = _element->GaussionPoint(j)[0];//_GaussPnt[j][0];
		        double eta = _element->GaussionPoint(j)[1];
		        double xi_ = _element->Global_x(xi, eta);
		        double eta_ = _element->Global_y(xi, eta);
                Dofs<2> temp({xi_,eta_});
		        a = a + _element->det_Jacobi(xi, eta) * _element->GaussionWeight(j) * _element->phi(xi, eta, i)* _func(temp.get_value());
            }
            _Rhs[_element->NdIdx(i)]+= a;
        }
    }
    //std::cout << _Rhs << std::endl;

}

void Poisson_2D::DealWithBoubdaryCondition()
{
    for(auto k:_mesh->Boundary())
    {
        Dofs<2> bnd_point = _mesh ->DofsofIndex(k);
	double bnd_value =_boundarycondition._bflist[0]->g()(bnd_point.get_value());
	 _Rhs[k] = bnd_value * _A.coeffRef(k,k);
	for(Eigen::SparseMatrix<double>::InnerIterator it(_A,k);it;++it)
	{
	    int row = it.row();
	    if(row == k)
		continue;
	    _Rhs[row] -= _A.coeffRef(k,row)* bnd_value;
	    _A.coeffRef(k,row) = 0.0;
	    _A.coeffRef(row,k) = 0.0;
	}
        
    }
    //std::cout << _A << std::endl;
    //std::cout << _Rhs << std::endl;
}

// void Poisson_2D::AssembleRightHandsTerm()
// {
//     int n_Dofs = _element->n_Dofs();
//     int n_GaussPnt = _element->n_GaussPnt();
//     auto _func  = _equation->RightHandsTermFunction();
//     for (int k = 0; k < _mesh -> n_element(); k++)
//     {
//         std::vector<int> Idx = _mesh->NodeofEle(k);
//         std::vector<Dofs<2>> temp(n_Dofs);
//         for(int i = 0; i < n_Dofs; i++)
//         {
// 	    //temp[i].GlobalIndex() = Idx[i];
// 	    temp[i] = _mesh ->DofsofIndex(Idx[i]);
//         }
//         _element->SetDofsList(temp);
//         for(int i = 1; i <= n_Dofs; i++)
//         {
//             double a = 0.0;
//             for (int j = 0; j < n_GaussPnt; j++)
//             {
// 		        double xi = _element->GaussionPoint(j)[0];//_GaussPnt[j][0];
// 		        double eta = _element->GaussionPoint(j)[1];
// 		        double xi_ = _element->Global_x(xi, eta);
// 		        double eta_ = _element->Global_y(xi, eta);
//                 Dofs<2> temp({xi_,eta_});
// 		        a = a + _element->det_Jacobi(xi, eta) * _element->GaussionWeight(j) * _element->phi(xi, eta, i)* _func(*temp);
//             }
//             _Rhs[_element->NdIdx(i)]+= a;
//         }
//     }
//     //std::cout << _Rhs << std::endl;
//     if(_boundarycondition.is_all_dirichlet())
//         return;
//     else
//     {
//         std::vector<BoundaryFunction<2>* > BFL = _boundarycondition.BFlist();
//         std::vector<std::vector<int> > BndEdge = _mesh->_BndEdge;
//         int num_of_bnd_edge = BndEdge.size();
//         int n_Dofs = _element ->n_Dofs();
//         int n_GaussPnt = _element ->n_GaussPnt();
//         for(int k = 0;k < num_of_bnd_edge;k++)
//         {
//             std::function<double(double*)> _func;
//             int bndtype = BndEdge[k][0];
//             if(bndtype == -2)
//             {
//             int _element_index = BndEdge[k][1];
//             std::vector<int> Idx = _mesh->NodeofEle(BndEdge[k][1]);
//             std::vector<Dofs<2>> temp(n_Dofs);
//             for(int i = 0; i < n_Dofs; i++)
// 	            temp[i] = _mesh ->DofsofIndex(Idx[i]);
//             _element->SetDofsList(temp);
//             int num1 = BndEdge[k][2];
//             int num2 = BndEdge[k][3];
//             Dofs<2> dof1 = _mesh->DofsofIndex(num1);
//             Dofs<2> dof2 = _mesh->DofsofIndex(num2);
//             bool is_horizontal = false;
//             bool is_vertical = false;
//             if(dof1[1] == dof2[1])
//                 is_horizontal = true;
//             else if(dof1[0] == dof2[0])
//                 is_vertical = true;
//             for(int i = 0;i < BFL.size();i++)
//             {
//                 BoundaryFunction<2> * t = BFL[i];
//                 for(int j = 0;j < t->Boundaryvector().size();j++)
//                 {
//                     Boundary<2> B = t->Boundaryvector()[j];
//                     if(B.is_on_boundary(dof1) && B.is_on_boundary(dof2))
//                     {
//                         _func = t->g();
//                         break;
//                     }
//                 }
//             }
//             if(is_horizontal)
//             {
//                 int c = dof1[1];
//                 bool is_bottom = (c == _mesh->RDomain()->LeftBottomPoint()[1]);
//                 for(int i = 1; i <= _element->n_Dofs(); i++)
//                 {
//                 double a = 0.0;
//                 for (int j = 0; j < n_GaussPnt; j++)
//                 {
// 		            double xi = _element->GaussionPoint(j)[0];
// 		            double eta = _element->GaussionPoint(j)[1];
// 		            double xi_ = _element->Global_x(xi, eta);
//                     Dofs<2> temp({xi_,double(c)});
//                     int c_ = is_bottom ? -1 : 1;
// 		            a = a + sqrt(_element->det_Jacobi(xi, eta)) * _element->GaussionWeight(j) * _element->phi(xi, double(c_), i)* _func(*temp);
//                 }
//                 _Rhs[_element->NdIdx(i)]+= a / 2.0;
//                 }
//             }
//             else if(is_vertical)
//             {
//                 int c = dof1[0];
//                 bool is_left = (c == _mesh->RDomain()->LeftBottomPoint()[0]);
//                 for(int i = 1; i <= _element->n_Dofs(); i++)
//                 {
//                 double a = 0.0;
//                 for (int j = 0; j < n_GaussPnt; j++)
//                 {
// 		            double xi = _element->GaussionPoint(j)[0];
// 		            double eta = _element->GaussionPoint(j)[1];
// 		            double eta_ = _element->Global_y(xi, eta);
//                     Dofs<2> temp({double(c),eta_});
//                     int c_ = is_left ? -1 : 1;
// 		            a = a + sqrt(_element->det_Jacobi(xi, eta)) * _element->GaussionWeight(j) * _element->phi(double(c_), eta, i)* _func(*temp);
//                 }
//                 _Rhs[_element->NdIdx(i)]+= a / 2.0;
//                 }
//             }
//             }
//         }
//     }
//     //std::cout << "neumann after" << std::endl;
//     //std::cout << _Rhs << std::endl;

// }

// void Poisson_2D::DealWithBoubdaryCondition()
// {
//     std::vector<BoundaryFunction<2>* > BFL = _boundarycondition.BFlist();
//     if(_boundarycondition.have_dirichlet())
//     for(auto k:_mesh->Boundary())
//     {
//         bool found_dirichlet_bnd = false;
//         Dofs<2> bnd_point = _mesh ->DofsofIndex(k);
//         double bnd_value ;
//         Boundaryinfo bctype;
//         Boundaryinfo bctype_for_n;
//         for(int i = 0; i < BFL.size();i++)
//         {
//             BoundaryFunction<2> * t = BFL[i];
//             if(!found_dirichlet_bnd)
//             {
//             bctype =BFL[i]->boundaryType();           
//             for(int j = 0;j < BFL[i]->Boundaryvector().size();j++)
//             {
//                 Boundary<2> B =  BFL[i]->Boundaryvector()[j];
//                 if(B.is_on_boundary(bnd_point))
//                 {
//                     if(bctype != 1)
//                     {
//                         bctype_for_n = bctype;
//                         break;
//                     }
//                     else
//                     {   
//                         bctype_for_n = bctype;
//                         bnd_value = BFL[i]->g()(*bnd_point);
//                         found_dirichlet_bnd = true;
//                     }
//                 }
//             }
//             }
//         }
//         if(bctype_for_n == 1)
//         {
//             _Rhs[k] = bnd_value * _A.coeffRef(k,k);
//             for(Eigen::SparseMatrix<double>::InnerIterator it(_A,k);it;++it)
//             {
//                 int row = it.row();
//                 if(row == k)
//                     continue;
//                 _Rhs[row] -= _A.coeffRef(k,row)* bnd_value;
//                 _A.coeffRef(k,row) = 0.0;
//                 _A.coeffRef(row,k) = 0.0;
//             }
//         }
//     }
//     //std::cout << _A << std::endl;
//     //std::cout << _Rhs << std::endl;
// }

void Poisson_2D::Solve()
{
    Eigen::ConjugateGradient<Eigen::SparseMatrix<double> > Solver_sparse;
    Eigen::VectorXd x1_sparse;
    Solver_sparse.setTolerance(1e-12);
    Solver_sparse.compute(_A);
    _Solution = Solver_sparse.solve(_Rhs);
    // std::cout << "solution is : " << std::endl;
    //std::cout << _Solution << std::endl;
}
double Poisson_2D::Error_H1Norm(std::function<double(double *)> real_solution,std::function<double(double *)> real_solution_x,std::function<double(double *)> real_solution_y)
{
    int n_Dofs = _element->n_Dofs();
    int n_GaussPnt = _element->n_GaussPnt();
    double h1error = 0.0,l2error=0.0;
    for(int k = 0;k < _mesh->n_element();k++)
    {
        double h1errorink = 0.0,l2errorink=0.0;
        std::vector<int> element_index = _mesh -> NodeofEle(k);
        std::vector<Dofs<2>> temp(n_Dofs);
        for(int i = 0; i < n_Dofs; i++)
	        temp[i] = _mesh ->DofsofIndex(element_index[i]);
        _element->SetDofsList(temp);
	for(int j = 0;j < n_GaussPnt;j++)
        {
            double u_j = 0.0;
            double u_h_j = 0.0;
            double xi = _element->GaussionPoint(j)[0];
            double eta = _element ->GaussionPoint(j)[1];
            double x = _element->Global_x(xi,eta);
            double y = _element ->Global_y(xi,eta);
            double guasspoint[2] = {x,y};
            u_j = real_solution(guasspoint);
            for(int i = 1;i <=n_Dofs;i++)
                u_h_j += _Solution(element_index[i - 1]) * _element ->phi(xi,eta,i);
            l2errorink += (u_j - u_h_j) * (u_j - u_h_j) * _element ->GaussionWeight(j) * _element ->det_Jacobi(xi,eta);
	}
        for(int j = 0;j < n_GaussPnt;j++)
        {
            double u_j_x = 0.0;
            double u_h_j_x = 0.0;
            double xi = _element->GaussionPoint(j)[0];
            double eta = _element ->GaussionPoint(j)[1];
            double x = _element->Global_x(xi,eta);
            double y = _element ->Global_y(xi,eta);
            double guasspoint[2] = {x,y};
            u_j_x = real_solution_x(guasspoint);
            for(int i = 1;i <=n_Dofs;i++)
                u_h_j_x += _Solution(element_index[i - 1]) * _element ->phi_x(xi,eta,i);
            h1errorink += (u_j_x - u_h_j_x) * (u_j_x - u_h_j_x) * _element ->GaussionWeight(j);
        }
        for(int j = 0;j < n_GaussPnt;j++)
        {
            double u_j_y = 0.0;
            double u_h_j_y = 0.0;
            double xi = _element->GaussionPoint(j)[0];
            double eta = _element ->GaussionPoint(j)[1];
            double x = _element->Global_x(xi,eta);
            double y = _element ->Global_y(xi,eta);
            double guasspoint[2] = {x,y};
            u_j_y = real_solution_y(guasspoint);
            for(int i = 1;i <=n_Dofs;i++)
                u_h_j_y += _Solution(element_index[i - 1]) * _element ->phi_y(xi,eta,i);
            h1errorink += (u_j_y - u_h_j_y) * (u_j_y - u_h_j_y) * _element ->GaussionWeight(j);
        }
        h1error += h1errorink;
    }
    h1error /= _mesh->n_dofs();
    l2error /= _mesh->n_dofs();
    return sqrt(h1error+l2error);
};
double Poisson_2D::Error_L2Norm(std::function<double(double *)> real_solution)
{
    int n_Dofs = _element->n_Dofs();
    int n_GaussPnt = _element->n_GaussPnt();
    double l2error = 0.0;
    for(int k = 0;k < _mesh->n_element();k++)
    {
        double l2errorink = 0.0;
        std::vector<int> element_index = _mesh -> NodeofEle(k);
        std::vector<Dofs<2>> temp(n_Dofs);
        for(int i = 0; i < n_Dofs; i++)
	        temp[i] = _mesh ->DofsofIndex(element_index[i]);
        _element->SetDofsList(temp);
        for(int j = 0;j < n_GaussPnt;j++)
        {
            double u_j = 0.0;
            double u_h_j = 0.0;
            double xi = _element->GaussionPoint(j)[0];
            double eta = _element ->GaussionPoint(j)[1];
            double x = _element->Global_x(xi,eta);
            double y = _element ->Global_y(xi,eta);
            double guasspoint[2] = {x,y};
            u_j = real_solution(guasspoint);
            for(int i = 1;i <=n_Dofs;i++)
                u_h_j += _Solution(element_index[i - 1]) * _element ->phi(xi,eta,i);
            l2errorink += (u_j - u_h_j) * (u_j - u_h_j) * _element ->GaussionWeight(j) * _element ->det_Jacobi(xi,eta);
        }
        l2error += l2errorink;
    }
    return sqrt(l2error);
};

SpMat Poisson_2D::A() const
{
    return _A;
}
Eigen::VectorXd Poisson_2D::Rhs() const
{
    return _Rhs;
}


Eigen::VectorXd Poisson_2D::solution()
{
    //std::cout<< "The solution of this possion equation is : "<<std::endl;
     return _Solution;
}

double Poisson_2D::LocalStiffMatrix(int i, int j)
{
    double a = 0;
    for (int k = 0; k < _element->n_GaussPnt(); k++)
    {
	double xi = _element->GaussionPoint(k)[0];
	double eta = _element->GaussionPoint(k)[1];
	double Innerproduct=0;
	for(int k = 0;k < _element->gradient(xi,eta,i).size();k++)
	{
	    Innerproduct += _element->gradient(xi,eta,i)[k]*_element->gradient(xi,eta,j)[k];
	}
	    a += _element->det_Jacobi(xi, eta) * _element->GaussionWeight(k) * Innerproduct;
    }
    return static_cast<double>(a);
};


#endif /* FEMSPACE_H */
